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Preface 


This  thesis  contains  approximate  solutions  to  the  lin¬ 
earized  equations  of  motion  for  the  spinup  of  an  ideal  dual¬ 
spin  spacecraft.  For  your  convenience,  a  list  of  the  sym¬ 
bols  used  in  the  development  of  these  approximations  can  be 
found  on  page  viii. 

I  would  like  to  thank  all  of  my  instructors  for  the 
clarity  and  understanding  with  which  they  taught  their  sub¬ 
jects.  I  would  especially  like  to  thank  two  faculty  mem¬ 
bers.  First,  I  am  grateful  to  Col  R.  L.  Bagley  for  making 
dynamics  really  understandable  so  that  it  became  my  favorite 
area  of  study.  Second,  I  thank  my  thesis  advisor,  Capt  C. 

D.  Hall,  for  his  timely  assignment  to  the  faculty  and  his 
interest  in  spacecraft  dynamics.  Without  him  I  would  have 
had  to  choose  a  thesis  topic  in  which  I  may  not  have  been 
interested. 


David  L .  Kinney 
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Abstract 


An  approximate  solution  is  obtained  for  the  linearized 
equations  of  motion  for  the  spinup  dynamics  of  a  dual-spin 
spacecraft.  The  approximation  is  obtained  from  the  WKB 
(Wentzel,  Kramers,  and  Brillouin)  method,  which  usually 
generates  a  divergent  series.  Comparison  between  the  WKB 
solution  and  the  numerically  integrated  solution  are  used  to 
study  the  effects  of  spinup  torque  magnitude  and  the  number 
of  WKB  terms  used  on  the  accuracy  of  the  approximation.  The 
WKB  method  generates  accurate  approximations  for  small  spin¬ 
up  torques  using  as  few  as  two  terms  in  the  series  when  no 
singularities  are  present  in  the  WKB  terms  during  the  spinup 
time  interval.  For  an  axisymmetric  spacecraft,  the  WKB  se¬ 
ries  truncates  after  one  term  and  provides  the  exact  solu¬ 
tion  to  the  equations  of  motion. 
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AN  APPROXIMATE  SOLUTION  FOR  THE  LINEARIZED  SPINUP 


DYNAMICS  OF  A  DUAL-SPIN  SPACECRAFT 


I .  Introduction 


Background 

Dual-spin  spacecraft  models  consist  of  a  platform  and  a 
rotor  initially  spinning  together.  To  get  the  platform  to 
despin,  angular  momentum  is  transferred  from  the  platform  to 
the  rotor  causing  the  rotor  to  spinup.  The  equations  of 
motion  for  these  dynamics  are  generally  only  solvable  by 
developing  an  approximate  solution  oi  oy  numerical  integra¬ 
tion.  The  approximate  solution  has  the  advantage  over  nu¬ 
merical  integration  in  that  if  the  state  of  a  system  needs 
to  be  known  at  a  particular  time  in  the  future,  only  that 
time  needs  to  be  substituted  into  the  approximation  to  ob¬ 
tain  the  answer.  The  numerical  method  must  be  integrated 
from  some  initial  time  when  the  state  is  known  to  the  time 
in  question.  In  1976,  Gebman  and  Mingori  published  such  an 
approximation  using  perturbation  methods  (reference  2) . 

Their  method  examined  flat  spin  recovery  which  will  not  be 
discussed  in  this  thesis. 

Problem  and  Scope 

This  thesis  contains  the  examination  of  the  dynamics  of 
a  dual-spin  spacecraft  using  a  perturbation  technique  called 
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the  WKB  (Wentzel,  Kramers,  and  Brillouin)  method.  The 
spacecraft  is  modeled  as  a  rigid  body  gyrostat  consisting  of 
a  rotor  and  a  platform  initially  spinning  together.  A  small 
constant  torque  is  applied  by  the  platform  on  the  rotor 
until  the  platform  is  no  longer  spinning  with  respect  to  an 
inertial  reference  frame.  The  equations  of  motion  for  this 
system  are  linearized  and  the  WKB  method  is  applied  to  ob¬ 
tain  an  approximate  solution.  Both  oblate  and  prolate  spin- 
up  are  examined.  The  linear  equations  of  motion  are  solved 
using  a  numerical  integrator.  The  WKB  approximation  and  the 
linear  solutions  are  compared  to  determine  the  accuracy  and 
limitations  of  the  WKB  solution.  Dimensionless  torques  of 
6=0.1,  0.01,  and  0.001  are  used  in  this  comparison  as  well 
as  from  one  to  six  terms  in  the  WKB  solution. 

Assumptions 

The  following  assumptions  are  made  in  this  problem: 

1 .  There  are  no  external  torques  on  the  spacecraft . 

2.  There  are  no  losses  due  to  friction. 

3.  Both  the  platform  and  the  rotor  are  rigid  bodies. 

4.  The  rotor  is  axisymmetric  about  its  spin  axis. 

5.  The  platform  is  asymmetric,  but  the  axisymmetric 

case  is  also  examined. 

6.  The  spinup  torque  is  constant  during  spinup  and 
zero  at  all  other  times. 
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General 


First  the  equations  of  motion  for  the  system  were  de¬ 
veloped.  The  form  of  these  equations  came  from  a  paper  by 
Hall  and  Rand  (reference  4)  and  its  related  PhD  dissertation 
by  Hall  (reference  3).  The  equations  were  linearized  and 
put  into  a  form  appropriate  for  application  of  the  WKB  meth¬ 
od.  The  WKB  approximation  was  then  developed  with  the  help 
of  a  mathematical  manipulation  program.  Next,  numerical 
integration  was  used  to  compare  the  linear  and  nonlinear 
solutions.  Then,  the  WKB  solution  was  compared  to  the  lin¬ 
ear  solution.  This  comparison  consisted  of  one  example 
spacecraft  for  oblate  spinup  and  one  for  prolate  spinup. 

The  spinup  torques  and  the  number  of  terms  in  the  WKB  ap¬ 
proximation  were  varied  during  each  comparison. 

Sequence  of  Presentation 

Chapter  II  contains  the  derivation  of  the  equations  of 
motion  of  the  dual-spin  spacecraft.  These  equations  are  put 
into  a  dimensionless  form,  linearized,  then  put  into  a  form 
that  can  be  used  for  the  WKB  method.  Spacecraft  spinup  is 
also  discussed.  Chapter  III  contains  a  description  of  the 
WKB  method  and  its  application  to  the  spinup  problem.  Ap¬ 
proximate  solutions  for  the  spinup  problem  are  developed. 
Chapter  IV  contains  the  comparisons  between  the  WKB  and 
linear  solutions  using  specific  examples.  Chapter  V  con¬ 
tains  the  major  conclusions  drawn  from  the  comparison  and 
recommendations  for  further  study. 


II. 


This  chapter  contains  the  description  of  the  spacecraft 
model  to  be  examined.  From  this  the  equations  of  motion  are 
derived,  put  into  terms  of  dimensionless  variables,  and  then 
linearized.  Finally,  spacecraft  spinup  and  the  significance 
of  the  moments  of  inertia  to  spinup  are  discussed.  For  a 
full  discussion  of  spacecraft  dynamics,  the  reader  is  re¬ 
ferred  to  the  book  Spacecraft  Attitude  Dynamics  by  Peter  C. 
Hughes  (Reference  5) . 

Spacecraft  Model  (3:24-26;4:641-642) 

The  dual-spin  spacecraft  to  be  examined  will  be  mod¬ 
elled  as  the  gyrostat  shown  in  Fig  1.  The  platform,  P,  and 
the  rotor,  R,  are  rigid  bodies  connected  by  a  rigid  shaft 
with  frictionless  bearings.  A  constant  spinup  torque  is 
applied  to  both  bodies  along  the  shaft.  There  are  no  exter¬ 
nal  torques  on  the  gyrostat  and  it  is  free  to  rotate  in 
space.  The  two  bodies  however  are  constrained  to  rotate 
about  the  axis  of  the  shaft  which  is  a  principal  axis,  @x, 
of  the  gyrostat.  To  keep  the  gyrostat  as  general  as  possi¬ 
ble,  the  platform  will  be  asymmetric  while  the  rotor  will  be 
axisymmetric  about  the  axis. 

Equations  of  Motion  (3 : 25-29 ; 4 : 642-643 ) 

In  Fig  1  the  unit  vectors,  £j,  j=l,2,3,  are  the  princi¬ 
pal  axes  of  the  gyrostat  and  are  fixed  in  the  platform.  The 
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vector,  h,  is  the  total  angular  momentum  of  the  gyrostat. 
Since  there  are  no  external  torques,  the  angular  momentum 
vector  is  fixed  in  the  inertial  reference  frame. 

Before  the  equations  of  motion  can  be  determined,  sev¬ 
eral  angular  velocities  must  be  defined.  The  angular  veloc¬ 
ities  of  the  platform  are  defined  as  0^,  o and  ©3  which 
correspond  to  the  three  principal  axes .  These  angular  ve¬ 
locities  are  relative  to  the  inertial  frame,  but  are  ex¬ 
pressed  in  terms  of  the  body-fixed  unit  vectors.  The  vari¬ 
able  (Og  is  the  angular  velocity  of  the  rotor  relative  to  the 
platform.  Therefore,  0^+GJg  is  the  angular  velocity  of  the 
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rotor  about  its  axis  of  symmetry  relative  to  the  inertial 
frame . 

Like  the  angular  velocity  of  the  platform,  the  angular 
momentum  of  the  gyrostat,  Ji,  can  be  written  in  terms  of  the 
body- fixed  coordinates  to  obtain  hx,  h2,  and  h3  along  the 
principal  axes .  Now  the  angular  momentum  can  be  written  as 
follows : 


hi 

Ji  o  o] 

©i 

.  N 

Itr 

ii 

*2 

*3. 

o  H"1 

H*  ° 

o  o 

_ i 

0>3_ 

+ 

1 - 

o  o 

1 _ 

(1) 


where  Ilf  I2,  and  I3  are  the  moments  of  inertia  about  the 
gyrostat's  principal  axes  and  I8  is  the  moment  of  inertia  of 
the  rotor  about  the  e3  axis.  I8  is  included  in  Ix  and  the 
transverse  moments  of  inertia  of  the  rotor  are  included  in 
I2  and  I2. 

Next,  the  derivative  of  Eq  (1)  is  taken  with  respect  to 
time,  t,  to  give 


dt 


dhx 

dt 

dh2 

’  0  -0)3  ©2  ‘ 

V 

dhx 

-©3*2  +©2*3 

dhj 

dt 

dh3 

.  dt  . 

+ 

©3  0-©! 

-©2  ©i  0 

^2 

*3. 

— -  +©3*.  -©,h, 
dt 

dh2 

-©2*3  +©1h2 

0 

0 

0 


(2) 


Note  that  Eq  (2)  is  set  equal  to  zero  since  there  are  no 
external  moments. 
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Two  more  terms  are  now  defined:  the  moment  of  inertia 

of  the  platform  about  is  Ip;  and  the  angular  momentum  of 
the  rotor  in  the  direction  is  ha  to  give 

I1^Ip  +  Ig  (3) 

and 


ha  =  Js(cos+o>1) 


(4) 


Combining  Eqs  (3)  and  (4)  and  the  relationship  for  hx  from 
Eg  (1)  gives 


®1* 


hx~ha 


(5) 


Then  using  Eq  (5)  and  the  relationships  for  h2  and  h3  from 
Eq  (1),  <&2,  and  can  be  eliminated  from  Eq  (2)  to  give 

the  equations  of  motion: 


dh3  _ 
dt 

dh2  = 
dt 

dh3  _ 
dt 

and  the  spinup  torque,  ga, 


is  defined  as 


(6) 

(7) 

(8) 


dha 

dt 


(9) 


A  first  integral  of  motion  may  be  obtained  by  noting 
that  there  are  no  external  moments  and  thus  the  angular 
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momentum  is  conserved.  This  yields  the  following  relation¬ 
ship: 

h2  -  hi  +  hi  +  h$  =  constant  ( 1° ) 


Dimensionless  Equations  (3 :30-31;4:643-644) 

The  equations  of  motion  can  be  further  simplified  by 
transforming  the  variables  into  dimensionless  form.  Using 
the  transformations  for  the  dimensionless  angular  momenta 

X,-^l  j- 1,2,3  (ID 

J  h 

the  dimensionless  moments  of  inertia 


ij  =  l-  *  j -1,2,3 


(12) 


the  dimensionless  angular  momentum  of  the  rotor  about  ej 


ha 


(13) 


the  dimensionless  time 


t  = 


h  t 


(14) 


and  the  dimensionless  torque  applied  by  the  platform  to  the 
rotor  about  §x 


e  - 

e'~P~ 


(15) 


Eqs  (6) -(9)  can  be  put  into  a  dimensionless  form: 
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•^1  “  ^2^3 


(16) 


X2-  (i3X1-|i)X3 


(17) 


i3=  (H-i2X1)X2 


(18) 


A  -e 


(19) 


where 


/•  )  «  d(  )  (20) 

1  }  ~zr 

Now  the  first  integral  of  motion  in  Eq  (10)  becomes 


Xl+Xl+Xl*  1  (2D 


Linearized  Equations 

Before  Eqs  (16) -(19)  can  be  linearized,  the  equilibrium 
points  must  first  be  found.  The  equilibrium  points  occur 
where  the  first  derivatives  with  respect  to  time,  t,  of  the 
momenta  are  all  equal  to  zero.  This  is  true  when  X2  and  X3 
are  equal  to  zero  for  the  case  e*0,  and  therefore  X3=±l  from 
Eq  (21)  (3:32-33).  Only  X3=+l  will  be  examined  since  this 
is  the  equilibrium  point  about  which  spacecraft  typically 
operate. 

Eqs  (16),  (17),  and  (18)  may  now  be  linearized  about 
the  equilibrium  point 
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where 


(22) 


(23) 


Substituting 


*1 

*2 

*3 


and 

X1  =  l  +xl 
X2  ~  x2 

*3=*3 


(24) 


(25) 


into  Eqs  (16),  (17),  and  (18)  and  noting  that  the  lower  case 
x's  are  small  perturbations  from  the  nominal  values  so  quad¬ 
ratic  terms  in  x  can  be  dropped  gives 


x3  =  0 

x2+  (U-i3)  x3  =  0 
x3+(i2-li)  x2  =  0 


(26) 

(27) 

(28) 
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For  an  axisymmetric  spacecraft  (i2=i3),  the  WKB  method  gives 
the  exact  solutions  for  Eqs  (26) -(28)  and  the  nonlinear 
equations.  This  will  be  shown  in  Chapter  IV. 

Next,  taking  the  derivatives  of  Eqs  (27)  and  (28)  gives 


x2  +  (n  -  i3 )  x3  +  |i  x3  =  0 


(29) 


x3  +  (i2-n)x2-|ix2  =  0 

Then  combining  Eqs  (27),  (28),  (29),  and  (30)  to  decouple 
the  variables  gives 


A  •  (3D 

x2  +_ - rx2  +  (|l-i3)  (Jl-i2)x2  =  0 

a3  "4 

A  .  (32) 

x3  + - rx3  +  ( i2 -jx)  (i3-n)x3  =  0 

22  “K 

In  order  to  use  the  WKB  method  to  solve  the  linear 
differential  equations  in  x2  and  x3,  Eq  (19),  its  first  in¬ 
tegral 

fi  =  et+H0  (33) 

and  the  change  of  variable 

x-et 

d()  _  d( ) 

~ar  ~3T  04) 

d2  (  )  .  c2  d2  (  ) 
dt^  dx2 
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must  be  substituted  into  Eq  (31)  to  give 


e 


2v" 

x2 


rr  %-rX2  +  ('C+^0-i3)  (X+^o-^)  x2  =  0 


(35) 


where 


"'’TFT 


(36) 


and  p-o  is  the  value  of  |i  at  t=0.  The  zero  subscript  nota¬ 
tion  on  any  other  variable  in  the  thesis  will  be  its  value 
at  some  initial  time.  The  only  exception  is  S0  which  will 
be  defined  later.  It  should  be  noted  that  Eq  (35)  is  the 
same  as  that  for  x3  if  the  subscripts  2  and  3  on  all  vari¬ 
ables  are  exchanged. 

Finally,  the  substitutions 


k2  -  |lo  ~  -^2 

k3  =  M-o  -  ^ 


(37) 


are  made  into  Eq  (35)  as  a  simplification  to  obtain 


2 

e2X2;  +  (x+k2)  (x+k3)x2  -  0  (38) 

x  +  k3 

Before  discussing  the  WKB  method,  a  few  words  need  to 
be  said  about  spinup  and  dimensionless  inertia  parameters. 


Spinup  (3:57-59) 

Spinup  of  an  ideal  dual-spin  spacecraft  (i.e.,  no  fric¬ 
tion)  usually  starts  when  the  platform  and  rotor  are  spin¬ 
ning  together  as  one  rigid  body  called  the  all -spun  condi¬ 
tion  where  cos=0.  The  spinup  motor  is  then  started  so  there 
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is  a  small  constant  torque  between  the  platform  and  rotor. 
Angular  momentum  is  transferred  from  the  platform  to  the 
rotor  such  that  the  rotor  is  spunup  and  the  platform  is 
despun.  The  motor  is  then  switched  off  when  (D^O.  Combin¬ 
ing  the  all-spun  condition,  (Os=0,  with  Eqs  (4)  and  (5)  gives 


(39) 


or  its  dimensionless  form 


(40) 


Combining  the  final  despun  condition,  0^=0,  and  Eq  (5)  gives 


ha  =  hi 


(41) 


or  its  dimensionless  form 


(42) 


The  linearized  counterpart  of  Eq  (40)  becomes 


(43) 


and  Eq  (42)  becomes 


H  =  1 


(44) 


Eqs  (43)  and  (44)  are  the  initial  and  final  spinup  condi¬ 
tions,  respectively,  for  linearized  spinup. 


Dimensionless  Moments  of  Inertia  (3 :33-34;4 : 644) 

According  to  Eq  (12),  the  relative  value  relationships 
between  the  actual  moments  of  inertia  and  their  dimension¬ 


less  counterparts  hold  true  such  that 


Ij>  !*<=>  ij->  ik  ,  j ,  k-1,  2,  3  (45) 

The  following  relationships  also  can  be  seen  from  Eq  (12): 


Ip>  Ij  <=>  ij  <  0  ,  j  -  2 ,  3 
Jp<  ij  >  0  ,  j  =  2,3 


(46) 


Now  the  dynamical  shape  of  the  gyrostat  can  be  defined. 

Note  that  in  Eqs  (6)-( 8)  only  i2  and  i3  appear.  So  the  pos¬ 
sible  dynamical  shape  of  the  gyrostat  depends  only  on  these 
two  parameters  and  oblate  and  prolate  are  now  defined  as 
follows : 

oblate:  i2<0  and  i3<0 

prolate:  i2>0  and  i3>0 

To  obtain  the  range  of  values  for  the  other  moment  of 
inertia,  i3,  Eqs  (3)  and  (12)  can  be  combined  to  give 


Is 

“3T 


i,  •  —  * 


'TP~rs 


(47) 


which  is  the  ratio  of  the  moment  of  inertia  of  the  rotor  to 
that  of  the  entire  gyrostat  about  the  axis .  This  means 
ix  can  take  on  the  values  of  0<i1<l  depending  on  the  rela¬ 
tionship  between  Ip  and  Is. 

Although  i3  does  not  appear  in  the  equations  of  motion, 
it  does  have  an  effect  on  the  spinup  dynamics.  In  the  pre¬ 
vious  section,  the  spin  up  was  seen  to  begin  where  the  ini¬ 
tial  value  of  was  equal  to  i3  for  the  linearized  case.  It 
also  limits  the  initial  rotor  spin  axis  to  the  major  or 
minor  axis  where: 
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major  axis:  i1>{i2,  i3} 


minor  axis:  i1<{i2,  i3}  . 

If  the  value  of  ix  falls  between  i2  and  i3,  the  all-spun 
condition  is  about  the  intermediate  axis  and  is  therefore 
unstable . 

Also,  the  physically  possible  values  of  i2  and  i3  are 
limited  by  i3.  The  inequalities 

^i<  ^2  +  ^3 

I2<J1+I3  (48) 

J3  <  Ix  +  I2 

must  hold  true  for  a  physically  possible  system.  For  more 
information  on  these  physical  limitations,  the  reader  is 
referred  to  reference  3  (3:176-180). 
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Ill .  WKB  Approximat ion 


This  chapter  contains  the  description  of  the  WKB  method 
and  how  it  is  applied  to  the  spinup  of  a  dual-spin  space¬ 
craft.  Approximate  solutions  are  obtained  for  oblate  spinup 
and  then  for  prolate  spinup. 

Background 

WKB  theory  is  a  method  for  obtaining  an  approximate 
solution  to  a  linear  differential  equation  where  the  highest 
derivative  is  multiplied  by  a  small  parameter  (1:484) .  It 
is  sometimes  referred  to  as  WKBJ  for  Wentzel,  Kramers, 
Brillouin,  and  Jeffreys  who  independently  used  this  method 
in  quantum  mechanics  during  the  1920s.  It  was  also  indepen¬ 
dently  discovered  by  many  others  so  that  it  is  also  some¬ 
times  called  Liouville-Green  or  phase  integral  (6:244)  . 

WKB  theory  contains  boundary-layer  theory  as  a  special 
case.  Boundary- layer  theory  will  not  be  discussed  here.  The 
reader  is  referred  to  reference  1  for  a  discussion  of  bound¬ 
ary-layer  theory.  The  WKB  approximation  starts  by  assuming 
solutions  of  the  form: 

x(T)  ~A(X)  exp[£^I]  ,  8-»0+  (49) 

where  S(x)  is  the  phase,  8  is  the  boundary  layer  thickness, 
and  A  (x )  is  an  amplitude  function.  The  symbol  is  used 
here  to  denote  the  asymptotic  nature  of  the  approximation. 

In  this  case,  asymptotic  means  as  8  gets  closer  to  0+  the 


difference  between  the  approximation  and  x(x)  becomes  much 
less  than  x(x) .  Eg  (49)  is  not  in  a  suitable  form  to  obtain 
an  asymptotic  approximation  because  there  is  an  implicit 
dependence  on  8  in  both  the  phase  and  amplitude  functions. 
This  implicit  dependence  can  be  made  explicit  by  expanding 
A ( X )  and  S(x)  as  a  power  series  in  8.  The  two  series  can 
then  be  combined  to  form  a  single  power  series  so  that  Eq 
(49)  becomes 

x(x)  -exp  8 Q  Sq(X) 

_  8  qrn  0 

from  which  an  approximate  solution  can  now  be  derived 
(1:484-486) . 

There  are  several  limitations  to  the  WKB  method.  The 
differential  equation  must  be  linear  and  its  highest  deriva¬ 
tive  must  be  multiplied  by  a  small  parameter  (1:484).  This 
is  the  reason  the  equations  of  motion  for  the  dual-spin 
spacecraft  were  put  into  the  form  of  Eq  (38) .  Also,  the  WKB 
series  within  the  exponential  of  Eq  (50)  is  not  convergent 
(1:493).  This  means  that  some  finite  number  of  terms  will 
give  the  best  approximation. 

WKB  Method 

This  section  contains  the  derivation  of  the  first  six 
terms  in  the  WKB  series  so  that  solutions  may  be  approximat¬ 
ed  for  x2  for  both  oblate  and  prolate  spinup.  The  steps  for 
this  method  are  taken  from  reference  1  (1:486-490).  The 
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derivations  of  all  the  WKB  terms  were  obtained  using  the 
symbolic  manipulation  program  Mathematica  2.0  for  SPARC  by 
Wolfram  Research,  Inc. 

The  first  step  to  obtain  the  first  six  terms  of  the  WKB 
series  is  to  take  Eq  (50)  and  its  derivatives 


x2  -  exp 


E  s*-1*, 

<3*0 


(  \ 
5 

(  \ 
5 

II 

E  &9_l5* 

exp 

E 

! 

lg*°  J 

J 

x"  = 


\2 


E 


<3*0 


+  E  s^s" 

<3*0 


r  \ 

5 

exp 

E  8’‘1  s. 

J 

L5-0  J 

(51) 


(52) 


(53) 


and  substitute  them  into  Eq  (38) : 


f  > 

2 

E  s’-1 4 

+ 

J 

V 

E  &r's” 

<3*0 


> 

r  \ 

exp 

E  8«-‘sa 

/- 

l®*0  J 

X  +  k, 


r  \ 

E 59-1 4 

exp 

E  59-1 

v«-°  J 

l9*0  J 

+  (X+k2)  (X+k3)  exp 


E  5<?_1 5<5 


yq*o 


=  0 


(54) 


After  eliminating  the  exponentials  and  expanding  Eq  (54), 
the  dominant  terms  are  found  to  be 


(x  +  k2)  (x  +  k3) 

e25o2 

“F“ 


(55) 
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For  these  terms  to  be  of  the  same  order  of  magnitude,  8  must 
be  of  the  same  order  as  e.  For  simplicity  8=e  is  chosen. 

The  next  step  is  to  make  the  substitution  8=e  into  Eq 
(54)  and  group  terms  by  powers  of  e: 


e°[  (x  +k2)  ( x  +  k3 )  +  So2]  +  e1 


'  ttSl  +  2S'0sl  *  si1' 


+  e- 


+  e3  s2  +  2  s{  S2'  +2  s'o  si  +  s2"  j 

J4  si2  +  si  +  2  si  si  +  2  So  si  +  s" 

^  3 

2  S2  S3'  +  -li-  S4  +  2  si  si  +  2  So  s'  +  S^  I 
X  +/c3  ^ 


) 


+  o(e6)  =  0 


(56) 


Now  the  terms  for  each  power  of  e  are  individually  set  equal 
to  zero  in  order  to  obtain  six  equations,  one  for  each  of 
the  derivatives  of  the  WKB  series  terms: 


So 


si 


si 


si 


±  [ -  (x  +  k2)  (X  +  k3)]1/2 
(X+J^r1  S0'  -  si' 

2SJ 

(x  +  ^j-1  s(  -  si2  -  si ' 


^  Sq 


(X+k-,)’1  Si  -  2  Si  Si  -  Si' 
2So 


(57) 

(58) 

(59) 

(60) 
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(61) 


si 


(X+k3)'1S: 


0/2  0  /  /  // 
02  “  z  03  —  03 


5s  = 


2  56 

( X  +  k3 )  -1  54  -  2  52  53  -  25(54'  ~  S[' 

IF0 


(62) 


Note  that  there  are  two  solutions  for  50'  and  thus  two  solu¬ 
tions  to  the  second  order  differential  equation  as  would  be 
expected.  To  determine  these  solutions  the  sign  of  the  term 
(t+k2)  (x+k3)  must  be  determined.  This  is  done  in  the  fol¬ 
lowing  sections  for  oblate  and  prolate  spinup. 

WKB  Oblate  Solution 

For  oblate  spinup,  both  i2  and  i3  are  negative  and  p.0 
is  assumed  to  be  positive.  This  means  the  derivative  of  the 
zeroth  term  in  Eq  (57)  is  imaginary  and  the  derivatives  of 
the  first  six  terms  of  the  WKB  series  may  be  written  as 
follows : 


Sq  =  ±i[(x+/e2)  ( x  +  k3 ) ] 


1/2 


si  = 


k2-k3 

4  (x+k,)  ( x  +  k, } 


c,=  +  i.  (k2-k2)  (7k2  +  5k3+12x) 
2  ±Z  32 [(x  +  k2)  (x+k3)]v^ 


c' =  3  (k3-k2)  (7k2  +8k2k3  +  5k3  +22k2T  +  18k3X  +20X2) 
3  64[(X+Jc2)  (X+k3)]1 


(63) 


(64) 


(65) 


(66) 
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si  (k2-k3)[{x+k2)  (x+k3)]~11/2  (1463ki+2163k2k3 

+  1989k2k2  *  11 05k2  +  655 2k\x  +  8304k2k3x 

+  5304/c|  T  +  10704)c2X2  +  9456k3X2  +  6720X3)  (67) 


Si  =t^(*2 -k3)[(X+k2)  (X  +Jc3)]"7  (707.k2  +1246)c|Jc3 

+  1392k|^3 +1130^2k33 +  565^3  +4074^|t 
+  6522Jcfk3X  +  6174Jc2Jcfx  +  3390Jc3x  +  9372)c2X2 
+  12696/c2/c3x2  +  8172/c32X2+10480Jc2X3  +  9680*:3X3  +  5040X4)  (68) 

Now  Eqs  (63) -(68)  can  be  integrated  to  obtain  the  first  six 
terms  of  the  WKB  series  for  oblate  spinup: 


=  ±i|  ^(X+Jc2)  (X+Jc3)]1/2 


T“  7 

-  ln{*:2  +  )c3  +  2x+2[(  x+&2  )  (  )]1/2} 


(69) 


5i  = 


In  (  x  +  )c3 )  -  In  (  x  +k2 ) 

- 3 - 


(70) 


S2  =  ^  (Ikl  -27k2k3  +  9k2k\  - 5k2  -  6k2x~36k2k3x  -  6Jt32X 

-24)c2x2 -24)c3x2 -16x3)  (k3-k2)'2[(x*k2)  (X+k3)]~3/2 


(71) 


(k2-k3)  Hk2+5k3+12x) 
64[(x  +  k2)  (X+k3)f 


(72) 
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Sa-±i 


[(x  +  fc2)  (x  +Jc3)]1/2  . 


1105 


9216  (Jc2 -/c3)  2  (%+k2) 


9216  (k2  -Jc3) 3  (x  +k2) 4 

15360  ()c2  -  )c3) 4  (X  +  )c2) 3 

1633 

127 

46080  U2-Jc3)b(X+/c2) 

2  5760  ()c2  -Jc3)  b  (X  +  )c2) 

1463 

*  2653 

9216  0c2-;c3)2(X+/c3)!j 

9216  (k2-k3)3{X+k3)i 

4973 

1  C  O  C  A  /  7*  7.  \  4  /  i  7.  \  3 

*  15767 

a  c  nan  t  \  \  ^ 

4223 


5760  (Jfe2-Jfe3)  *»  (X 


+  *3>  ] 


(73) 


Sc  = _ _ _  (707)c2  +  1071)e2Jc3  +  1017)c2)e3 

2048[(X  +  Jk2)  (X  +  Jc3)]fc  3  2 

+  565Jc3  +3192)c|x  +4176Jt2Jc3X  +2712)c|x  +  5280)e2x2 

+  4800)c3X2+3360x3)  (74) 

Since  only  the  even  numbered  terms  are  imaginary  for 
the  oblate  spinup,  the  WKB  series  can  be  divided  into  imagi¬ 
nary  and  real  parts: 

E  e'-'s.-E  (75> 

g»0  ««0  n«l 


where 

Q  =  M  +  N 
M  £  N  £  M+l 

The  m  summation  contains  the  imaginary  terms,  while  the  n 
summation  contains  the  real  terms.  Using  Eq  (75)  and  noting 
that  the  two  WKB  series  solutions  are  complex  conjugates  of 


22 


each  other,  the  general  WKB  solution  for  oblate  spinup  may 
be  written  as 


x2  =  q  exp  52  e2”1'1  s2a  +  e2  (n'1)  s2n.x 

mm0  n* 1 


+  C2exp  -52  e2fl"1S2fll  +  52  e2(/I"1,  S2n.x  {76 

jn-0  n-1  . 


where 


Cx  =  constant 
C2  =  constant 

The  solution  for  x3  could  just  as  easily  be  written  by  not¬ 
ing  the  symmetry  in  Eq  (35)  that  was  mentioned  in  Chapter 
II .  Cx  and  C2  can  be  obtained  by  imposing  the  following 
initial  conditions  on  Eq  (76) : 

x2lx0)  =x20 

(77 

x2  (T0)  =  X20 


The  result  is 


x20+x20  £  ^2m  1  &2  m.  0-T  ^  ^2n-l,0 

.m-0  n-1 


M 

c2m-l  qf 
e  *2m,  0 

M 

~  —  -  . . .  ■ 

N 

E 

exp 

E 

p2m-X  o  +  \  ■*  e2(n-l)  o 

fc  *2m,  0  *  2Lt  fc  *211-1,0 

m-0 

) 

^i-0 

n-1  J 

x20~x2oT^  e2fl,'1S2ln,0  +  ^  e2*13'11  S2n.1(0 


Ef2(77-1)  a 

e  *2n-l,0 


2m-l  c*' 


where  the  notation  Sq0  means  the  initial  value  of  Sq.  Using 
the  results  of  Eqs  (76),  (78),  and  (79),  and  eliminating  the 
imaginary  numbers  by  substituting  the  identities 


ft  _  exp  (i8)  +  exp  ( -i0) 

=  exp  (i0)  -  exp  (-20) 
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(80) 


gives 


x2  “exp 


X)  e2(n  15  (s2n_i  -  s2n_1(0) 

77  “1 


r  -f 

M 

x20  COS 

-i£  e2-1(s2„-s2„,0) 

. 

in*0 

N 

^  r 

*20-*ao£  e2,n”1)5,2/n-i.( 


n"l 


M 


-i  e2”"1  s2in  o 

m»0 


sin 


“i]£  e2m  1(S2m  ~  S2m,o) 

ni”  0 


(81) 


which  is  the  WKB  solution  for  the  oblate  spinup  problem. 


WKB  Prolate  Solution 

For  prolate  spinup,  both  i2  and  i3  are  positive  and  the 
magnitude  of  n0  is  assumed  to  be  less  than  the  magnitude  of 
either  i2  or  i3.  So  from  Eg  (37)  both  k2  and  k3  are  nega¬ 
tive.  This  means  the  derivative  of  the  zeroth  term  in  Eq 
(57)  can  be  either  real  or  imaginary.  The  prolate  spinup 
can  therefore  be  divided  into  three  regions,  each  with  a  WKB 
approximation.  To  determine  these  regions,  the  assumption 
k2<k3  is  made  so  the  three  regions  are  as  follows: 
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Region  1 :  0^X<-k3 

Region  2 :  -k3<X<-k2 

Region  3 :  -k2<x£l-|i0 

It  should  be  noted  that  Eqs  (64) -(68)  have  singular  points 
at  X=-k2  and  x=-k3  so  that  solutions  near  these  points  may 
not  be  valid.  No  generality  is  lost  when  making  the  assump¬ 
tion  k2<k2  since  if  k2>k3  then  a  coordinate  transformation 
would  bring  the  situation  back  to  the  original  assumption. 

Region  1  Solution .  Since  x<-k2  and  x<-k3  in  region  1, 
the  term  (x+k2)  (X+k3)  is  always  positive.  Therefore,  the 
derivative  of  the  zeroth  WKB  term  is  imaginary.  This  yields 
the  same  solution  as  for  oblate  spinup. 

Region  2.  Solution.  Since  x>-k3  and  x<-k2  in  region  2, 
the  term  (X+k2)  (X+k3)  is  always  negative.  Therefore,  the 
derivative  of  the  zeroth  WKB  term  is  real  and  the  deriva¬ 
tives  of  the  first  six  terms  of  the  WKB  series  may  be  writ¬ 
ten  as  follows: 

Sq  =±[-(x+k2)  (x+k3)]1/2  (82) 

S>1  "tTxT 4iT^-7c3T  (83> 


=± 


(k2-k3)  (7k2  +  5k3  +  12x) 
32[- (X  +k2)  (x+*3)jiya 


(84) 


si 


3  (k3-k2)  (7k2  +8k2k3  +  5k3  +  22k2x  +  18k3X  +20x2) 
64[(x  +k2)  (X+k3)]4 


(85) 
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si  (/c2-.fc3)[-(X  +  ;c2)  (X  +  /c3)]"U/2  (1463^1  +  2163^3 

+  1989k2k3  +1105^c3  +  6552Jc2X  +  8304Jc2)c3x 

+  5304/c3X  +  10704)c2X2  +  9456Jc3X2  +  6720X3)  (86) 

Si  =T^(k2-k3)[(X+k2)  (X  +  )c3)]'7  (707Jc2  +  1246k2k3 

+ 13  92Jc|jt|  +1130k2k33  +  565^3  +  4074*:23X 
+  6522k2k3x  +  6114k2k3X  +  3390k3X  +  9372Jc|x2 
+  12696Jt2Jc3X2  +  8172.k|x2  +  10480.k2X3  +  9680Jc3X3  +  5040x4)  (87) 

Now  Eq  (82) -(87)  can  be  integrated  to  obtain  the  first  six 
terms  of  the  WKB  series  for  prolate  spinup  in  region  2 : 


;0  =  ±{f*L^4V(t*k2)  <x»*3>]*/2 


( k3  ~k2  ) 2  .  ( k2+k3+2x 

5  \k2-k3\  Jj 


(88) 


S,  = 


In  (  x  +k3 )  -  In  ( x  +k7 ) 


(89) 


S2  =  ±  (Ikl  -21k\k3  +  9k2k3  -  5 k]  -  6 kjx  - 3 6k2k3x  -  6k%X 

-24/c2X2-24/c3x2-16x3)  (k3-k2)  '2[- ( x  +k2)  (x+/c3)]"3/2  (9Q) 


S->  - 


_  (k2-k3)  {lk2  +5k3  +12x) 


64[(X  +k2)  (X  +  )c3)]3 


(91) 
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SA=±[-  (X+k2)  (x  +  k3)f'2  . 


1547 
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5c  = _ -3— k2 _ -  (707/c|  + 1071^1  Ar3  +  1017^2Jt32 

5  2048[(X+*2)  (X+*3)f  3  2 

+  565Jc3  +  3192.kf  X  +  All  6k2k3X  +  2712.kf  X  +  5280Jc2X2 

+4800*3X2+3360x3)  (93) 

As  with  the  oblate  spinup  problem,  the  WKB  series  for 
region  2  of  prolate  spinup  can  be  divided  into  even  and  odd 
terms  except  all  terms  are  real  this  time: 

M+N  M  N 

E  e<7'ls?=E  +  E  e2(n-1>52l,_1  (54) 

orO  m«0  n»i 

Using  Eg  (94)  and  noting  that  the  two  WKB  series  solutions 
only  differ  by  a  factor  of  negative  one  in  their  even  terms, 
the  general  WKB  solution  for  prolate  spin-up  in  region  2  may 
be  written  as 
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x2  =  C,  exp 


fli*0 


E  e-2a~1  S2w  +  E  e2(n'1)52n-l 

n«l 


+  C2exp 


"E  +  E  e2(n*1>  s2j3_a 

m«0  n«l 


(95) 


Cx  and  C2  can  be  obtained  by  imposing  the  following  initial 
conditions  on  Eg  (95) : 


x2(x0)  =x20 

X2'  <T0)  =X20 


(96) 


The  result  is 
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(98) 


Using  the  results  of  Egs  (95),  (97),  and  (98)  and  the  iden¬ 
tities 


sinh  (8)  >  e*P(8>  -«*P<-»> 

Z  (99) 

cosh (6)  =  exp(6)  +  exp(-0) 


gives 
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x2  =  exp 
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(100) 


which  is  the  WKB  solution  for  region  2  of  the  prolate  spinup 
problem.  Note  how  this  solution  parallels  the  oblate  solu¬ 
tion  in  Eq  (81) . 

Region  3.  Solution.  Since  X>-k2  and  x>-k3  in  region  3, 
the  term  (X+k2)  (x+k3)  is  always  positive.  Therefore,  the 
derivative  of  the  zeroth  WKB  term  is  imaginary.  Like  the 
solution  for  region  1,  this  yields  the  same  solution  as  for 
oblate  spinup. 

Now  that  the  WKB  approximations  have  been  derived,  they 
can  be  compared  to  the  linear  solutions  to  determine  their 
accuracy.  This  is  done  in  the  next  chapter. 
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This  chapter  contains  the  results  of  the  study  of  the 


WKB  method.  After  a  brief  description  of  the  procedure 
used;  a  comparison  of  the  nonlinear  and  linear  solutions  for 
oblate  spinup  is  presented.  Then,  the  WKB  solution  is  com¬ 
pared  to  the  linear  solution  for  oblate  spinup.  Next  is  a 
presentation  of  the  comparison  of  the  nonlinear  and  linear 
solutions  for  prolate  spinup.  Finally,  a  comparison  of  the 
WKB  and  linear  solutions  for  prolate  spinup  is  presented. 

Procedure 

All  of  the  numerical  integrations  of  the  nonlinear 
equations  of  motion,  Eqs  (16) -(19),  and  linear  equations  of 
motion,  Eqs  (26) -(28),  were  accomplished  using  a  fourth 
order  predictor-corrector  algorithm  developed  by  the  numeri¬ 
cal  analyst  Hamming.  This  algorithm  was  used  in  a  computer 
program  coded  in  FORTRAN  (7:107-113).  A  sample  of  the  type 
of  program  used  is  in  the  Appendix.  The  comparison  between 
different  solutions  was  done  in  the  same  program  and  the 
data  output  to  a  file.  This  file  was  then  imported  to  a 
spreadsheet  program  to  develop  the  graphs  in  this  chapter. 

Linear  vs  Nonlinear  Oblate  Solutions 

This  section  contains  an  examination  of  the  accuracy  of 
the  linear  solution  as  X1  gets  farther  from  the  equilibrium 
point.  The  effects  of  asymmetry  are  also  explained. 
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In  Fig  2,  the  magnitude  of  the  maximum  relative  error 
and  the  maximum  phase  error  between  the  nonlinear  and  linear 
solution  for  axisymmetric  (i2=i3)  oblate  spinup  is  plotted 
as  a  function  of  the  value  of  X3  for  several  values  of  e. 


Figure  2.  Linear  Error  in  x2  for  i2=i3=-0.5 


The  maximum  error  is  relative  to  the  maximum  amplitude  of 
x2,  a  constant  (i.e.,  the  magnitude  of  the  linear  solution 
minus  the  nonlinear  solution  all  divided  by  the  maximum 
value  of  x2) .  Spinup  is  taken  to  start  at  H=0  and  end  at 
ji=l  to  obtain  the  worst  case  error.  Fig  2  shows  that  for 
each  e  the  error  gets  worse  as  X3  departs  from  the  X3=l  val¬ 
ue  at  the  equilibrium  point.  Also,  smaller  e's  require  that 
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Xj  be  closer  to  the  equilibrium  point  because  as  e  decreases 
by  a  factor  of  ten,  the  spinup  time  increases  by  a  factor  of 
ten  which  gives  more  time  for  error  to  build  up. 

Fig  3  is  the  same  type  of  graph  as  Fig  2  except  only 
one  value  of  e  and  three  values  of  i2=i3  are  plotted.  Fig  3 


shows  that  the  smaller  the  inertia  parameter  the  larger  the 
error  for  a  given  Xx . 

Fig  4  shows  the  error  between  the  nonlinear  and  linear 
solutions  (i.e.,  linear  minus  nonlinear)  for  the  following 
conditions : 

Xx  =  0.9999 
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This  is  a  highly  asymmetric  body  and  departure  from  the 
nonlinear  solution  would  be  expected.  From  Fig  3  for 
X^O.9999  and  i2=i3=-0.9,  the  relative  error  is  approximate¬ 
ly  0.009.  Multiplying  the  relative  error  by  the  maximum 
value  of  x2  calculated  from  Eq  (21)  gives  an  error  of  ap¬ 
proximately  0.00013  for  the  axisymmetric  gyrostat.  The 
maximum  error  from  Fig  4  is  approximately  0.00023  which  is 
the  same  order  of  magnitude  as  for  the  axisymmetric  case, 
but  larger.  This  means  that  for  asymmetric  gyrostats,  the 
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value  of  Xj  must  be  closer  to  one  than  for  an  axisymmetric 
gyrostat  {with  moments  of  inertia  the  same  as  the  lesser  of 
i2  or  i3)  to  produce  the  same  linear  error.  Linear  error 
may  also  be  reduced  by  making  the  gyrostat  less  asymmetric. 
From  this  information,  the  range  of  Xx  can  be  determined  for 
a  particular  gyrostat  so  that  the  linear  solution  approxi¬ 
mates  the  nonlinear  solution. 

WKB  vs  Linear  Oblate  Solutions 

This  section  contains  the  results  of  a  comparison  of 
the  WKB  and  linear  solutions  for  oblate  spinup.  The  follow¬ 
ing  conditions  were  used  throughout  this  section  for  compar¬ 
ison: 

Xj*  0.9999 
H0=  0.1 
i2=-0 .6 
i3=-0 .4 

The  comparisons  were  made  using  three  values  for  e  (0.1, 
0.01,  and  0.001)  and  from  one  to  six  WKB  terms.  The  value 
of  Xx  and  the  choice  of  i2  and  i3  should  give  a  good  linear 
approximation  for  all  e .  The  program  used  to  run  this  com¬ 
parison  is  in  the  Appendix. 

Fig  5  shows  the  linear  solution  and  WKB  approximation 
for  e=0.1  and  one  WKB  term.  For  all  values  of  e,  the  ap¬ 
proximation  using  one  WKB  term  is  the  only  case  to  show  any 
significant  error.  Using  only  one  WKB  term  only  gives  an 
estimate  of  the  phase  while  the  amplitude  is  held  constant. 
The  time  varying  amplitude  correction  does  not  appear  until 
a  second  WKB  term  is  added  to  the  approximation.  The  third 
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Figure  5.  Linear  and  WKB  Solutions  for  e=0.1  and  One  WKB 


Term 

term  mainly  adjusts  the  varying  phase  again.  Then  the 
fourth  term  adjusts  the  varying  amplitude,  etc.  This  is 
evident  from  Eq  (81) . 

Fig  6  is  a  summary  of  all  of  the  coirparisons .  It  shows 
the  error  in  the  last  maximum  (or  minimum)  of  the  WKB  solu¬ 
tion  relative  to  the  last  maximum  (or  minimum)  of  the  linear 
solution  while  varying  the  number  of  WKB  terms.  Fig  6  only 
contains  only  amplitude  error.  It  shows  that  there  is  lit¬ 
tle  increase  in  accuracy  of  the  amplitude  after  two  terms 
for  e=0 .1  and  0.001.  But  there  is  still  significant  in¬ 
crease  in  accuracy  for  e=0.01  through  the  sixth  WKB  term. 
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Figure  6.  Relative  Error  in  the  Last  Maximum  of  x2  as  e  and 
the  Number  of  WKB  Terms  Change 


Fig  7  shows  the  maximum  error  in  the  WKB  solution  over 
the  entire  spinup  relative  to  the  maximum  magnitude  of  the 
linear  solution  during  the  spinup.  This  error  includes  the 
effects  of  both  magnitude  and  phase.  It  indicates  that  for 
e=0.1,  the  WKB  solution  eventually  diverges  from  the  linear 
solution  or  a  least  oscillates  around  some  final  value  as 
the  number  of  terms  increases.  For  e=0.01,  there  is  a 
steady  improvement  in  the  WKB  approximation  until  after  the 
fifth  term.  To  tell  whether  it  will  reach  some  final  value 
or  not  would  require  the  evaluation  of  more  WKB  terms.  For 
e=0.001,  the  error  of  the  WKB  solution  reaches  a  constant 
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Figure  7 .  Maximum  Relative  Error  in  x2  as  e  and  the  Number 
of  WKB  Terms  Change 

value  after  the  fourth  term.  The  apparent  oscillatory  be¬ 
havior  of  the  e=0.1  error  is  probably  due  to  e  being  not 
much  less  than  one  as  WKB  theory  requires.  The  higher  order 
terms  are  multiplied  by  increasing  powers  of  e.  To  reach  a 
final  value  quickly,  e  must  be  much  less  than  one.  This  is 
why  the  approximation  with  e=0.001  reaches  its  final  value 
so  quickly.  After  three  terms,  the  maximum  relative  error 
for  e=0.01  is  less  than  for  e=0.001.  However,  if  the  spinup 
for  e=0.001  is  examined  over  the  same  period  of  time  as  the 
total  spinup  for  e=0.01,  error  in  the  £=0.001  approximation 
is  less. 


This  example  demonstrates  that  the  WKB  solution  is  very 
good  approximation  to  linear  oblate  spinup  dynamics .  Except 
for  the  e=0.1  case,  relative  errors  of  less  than  1%  can  be 
obtained  with  just  two  terms. 

The  effect  of  changing  the  inertia  parameters  was 
briefly  examined.  The  cases  investigated  had  the  same  con¬ 
ditions  as  the  previous  oblate  example  except  the  inertia 
parameters  were  varied.  Only  £=0.01  and  two  and  four  terms 
were  used.  Table  1  shows  the  results.  It  lists  the  inertia 
parameters  used  and  the  approximate  factor  by  which  the 
error  increased  over  the  i2=-0.6,  i3=-0.4  case. 


Table  1 


Error  Increase  as 

Inertia  Parameters 

are  Varied 

la 

i3 

2  TERMS 

4  TERMS 

-0.9 

-0.4 

1 

1 

-0.5 

-0.4 

1 

1 

-0.6 

-0.5 

0.5 

0.5 

-0.6 

-0.2 

5 

70 

-0.9 

-0.1 

10 

100 

-0.6 

-0.1 

10 

100 

-0.2 

-0.1 

10 

1000 

The  trend  seems  to  indicate  that  inertia  parameters  that  are 
closer  to  zero  than  those  used  in  the  baseline  case  will 
produce  significantly  larger  errors  than  the  baseline  case. 
This  is  most  likely  due  to  spinup  starting  closer  to  the 
singularity  at  c=i3-n0=-k3  in  the  WKB  solution. 
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Nonlinear  vs  Linear  Prolate  Solution 

This  section  contains  a  brief  look  at  two  of  the  fac¬ 
tors  that  cause  the  linear  solution  to  diverge  from  the 
nonlinear  solution  for  prolate  spinup. 

Fig  8  shows  a  comparison  of  the  nonlinear  and  linear 
solutions  for  the  following  conditions: 

e  =0.01 
X10=  0.9999 
i2  =  0.6 

ij  =  0.4 
=  0-1 


Figure  8.  Nonlinear  and  Linear  Solutions  for  e=0.01 

The  graph  has  been  divided  into  the  three  regions  defined  in 
Chapter  III.  The  above  conditions  cause  the  linear  solution 
to  deviate  from  the  nonlinear  solution.  There  is  an 
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overshoot  of  the  first  maximum  in  region  3  which  causes 
errors  in  both  the  extreme  amplitudes  and  in  the  phase  of 
the  oscillation  through  the  remainder  of  the  spinup. 


Fig  9  demonstrates  the  effect  of  changing  only  the 
initial  value  of  X1  from  the  conditions  of  Fig  8.  The  value 


of  X10  is  now  0.999  which  is  farther  from  the  equilibrium 
point.  The  overshoot  of  the  first  maximum  in  region  3  has 
now  increased.  This  causes  more  error  in  the  extreme  ampli¬ 
tude  value  and  the  phase  over  that  seen  in  Fig  8. 

Fig  10  demonstrates  the  effect  of  making  the  spacecraft 
more  asymmetric  than  the  case  in  Fig  8 .  The  value  of  i3  has 
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Figure  10.  Nonlinear  and  Linear  Solutions  for  i3=0.3 

now  been  changed  to  0.3  which  causes  a  deviation  from  the 
nonlinear  solution  which  is  even  greater  than  the  case  in 
Fig  9. 

Using  figures  similar  to  Figs  8-10,  it  can  be  shown 
that  as  e  increases  the  value  of  X10  must  be  closer  to  the 
equilibrium  point  and  the  spacecraft  must  be  less  asymmetric 
for  the  linear  solution  to  better  approximate  the  nonlinear 
solution. 

It  appears  that  region  2  is  where  exponential  growth 
takes  place  in  the  linear  solution  which  agrees  with  the 
approximate  solution  in  Eq  (100) .  Any  error  in  the  linear 
solution  will  rapidly  grow  in  region  2  and  will  become  pro- 
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nounced  if  region  2  is  too  wide.  As  e  increases,  the  time 
spent  in  region  2  increases,  which  will  all  amplify  the 
error . 

WKB  vs  Linear  Prolate  Solutions 

This  section  contains  the  results  of  a  comparison  of 
the  WKB  and  linear  solutions  for  prolate  spinup.  The  fol¬ 
lowing  conditions  are  used  throughout  this  comparison: 

Xx=  0.99999 
^0=  0.1 
i2=  0.6 
i,=  0.5 

Comparisons  are  done  using  three  values  for  £  (0.1,  0.01, 
and  0.001)  and  from  one  to  six  WKB  terms.  Although  the 
conditions  are  just  out  of  the  linear  region  for  e=0.001, 
the  comparison  still  shows  of  how  well  the  WKB  method  ap¬ 
proximates  the  linear  solution.  The  actual  comparisons  are 
divided  into  the  three  regions  for  prolate  spinup  defined  in 
Chapter  III. 

Region  .1  Comparison.  Region  1  of  prolate  spinup  starts 
at  T=0  and  ends  near  t=0.4  which  is  a  singularity  for  all 
but  the  zeroth  term  of  the  WKB  series  (S0)  .  Therefore,  one 
comparison  between  WKB  and  linear  solutions  could  be  as 
shown  in  Fig  11.  Fig  11  is  a  plot  that  shows  how  close  the 
WKB  solutions  get  to  the  t=0.4  point  before  they  deviate 
from  the  linear  solution  by  1%  relative  to  x20.  Fig  11 
shows  that  for  e=0.1  the  WKB  solution  with  one  term  came 
closest  to  the  singularity.  For  e=0.01,  four  terms  came 
closest.  For  e=0.001,  there  was  no  substantial  change  after 
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Figure  11.  WKB  vs  Linear  Comparison  for  1%  Relative  Error 


Criterion  in  Region  1 

the  third  terms,  however  the  solution  with  four  terms  was 
slightly  better  than  the  others.  Using  this  criterion, 
small  e's  were  better  than  larger  ones. 

Another  comparison  is  shown  in  Fig  12.  It  plots  the 
maximum  error  between  the  WKB  and  linear  solutions  relative 
to  x20  as  a  function  of  the  number  of  WKB  terms  used  for  the 
three  e's.  This  error  was  only  measured  between  T=0  and 
1-0.2  to  try  to  eliminate  the  error  caused  by  the  singulari¬ 
ty  at  1=0.4.  Fig  12  shows  that  on  this  interval,  the  best 
approximation  for  e=0.1  is  still  one  WKB  term.  For  e=0.01, 
five  terms  is  now  best.  For  e=0.001,  six  terms  shows  the 
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Figure  12.  Relative  WKB  Error  in  First  Half  Region  1 
least  error. 

Fig  13,  14,  and  15  are  three  examples  that  show  the 
comparison  of  the  WKB  and  linear  solutions  for  prolate  spin- 
up  in  region  1.  Fig  13  is  the  comparison  for  e=0.1  for  one 
WKB  term  which  was  the  best  approximation  for  both  the  Fig 
11  and  12  criteria.  Fig  14  is  the  comparison  for  £=0.01  and 
four  WKB  terms  which  was  the  best  for  the  Fig  11  criterion. 
And  finally.  Fig  15  is  the  comparison  for  £=0.001  and  four 
terms  whiun  was  also  the  best  for  the  Fig  11  criterion. 

The  two  comparison  criteria  in  Figs  11  and  12  were 
completely  arbitrary  and  other  similar  comparison  criteria 


could  have  been  used,  but  overall  the  smaller  the  e  used  the 
better  the  approximation  of  the  linear  solution. 
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Figure  15.  WKB  and  Linear  Solutions  in  Region  1  for  e=0.001 
and  Four  Terms 


Region  2_  Comparison.  Region  2  comparison  is  more  com¬ 
plicated  than  for  region  1.  Region  2  is  bounded  by  two 
singularities  in  the  WKB  solution  at  x=0.4  and  T=0.5.  Also, 
an  initial  condition  for  the  WKB  solution  in  region  2  can 
not  be  obtained  from  the  WKB  solution  in  region  1 .  The 
initial  condition  used  for  the  WKB  solution  in  region  2  came 
from  numerically  integrating  the  linear  equations  of  motion 
out  to  the  point  x=0.45  which  is  the  half  way  point  between 
the  singularities.  This  point  was  chosen  to  try  to  minimize 
the  effects  of  the  singularities  on  the  WKB  approximation  in 
this  region.  The  WKB  solutions  were  then  propagated  back¬ 
ward  and  forward  from  the  x=0.45  point  and  compared  to  the 
linear  solution. 

Fig  16  shows  the  same  type  of  graph  as  Fig  11.  Fig  16 
is  a  plot  that  shows  how  close  the  WKB  solutions  get  to  the 
x=0 . 4  and  0.5  points  before  they  deviate  from  the  linear 
solution  by  1%  relative  to  the  magnitude  of  the  extreme 
value  of  the  linear  solution  in  region  2.  Since  the  extreme 
values  are  different  for  each  e,  care  should  be  taken  when 
trying  to  compare  one  e  to  another.  The  closer  the  points 
are  to  the  top  and  bottom  of  the  graph  ,  the  better  the 
approximation  is.  For  e=0.1,  the  WKB  solution  with  one  term 
was  the  best.  For  e=0.01,  two  WKB  terms  were  best.  Two 
terms  were  also  best  for  6=0.001  case  based  on  this  criteri¬ 
on. 

Figs  17,  18,  and  19  are  examples  of  the  comparisons 
between  the  WKB  and  linear  solutions  for  prolate  spinup  in 
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Figure  16.  WKB  vs  Linear  Comparison  for  1%  Relative  Error 
Criterion  in  Region  2 


region  2.  Fig  17  is  the  best  approximation  for  6=0.1  (one 
WKB  term).  Fig  18  is  the  best  approximation  for  e=0.01  (two 
WKB  terms).  Fig  19  is  the  best  approximation  for  6=0.001 
(two  WKB  terms)  . 

The  WKB  solutions  in  region  2  are  not  as  accurate  as 
those  for  region  1.  If  region  2  were  larger,  the  WKB  solu¬ 
tion  would  probably  have  been  better  than  it  was,  but  this 
would  have  made  the  linear  solution  deviate  farther  from  the 
nonlinear  solution. 
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Figure  17.  WKB  and  Linear  Solutions  in  Region  2  for  e=0.1 
and  One  Term 
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Figure  19.  WKB  and  Linear  Solutions  in  Region  2  for  e=0.001 
and  Two  Terms 


Region  3.  Comparison .  In  region  3  the  initial  condi¬ 
tions  for  the  WKB  approximation  can  not  be  obtained  from 
region  2  because  of  the  singularity  at  x=0.5.  So  as  in 
region  2,  the  linearized  equations  of  motion  were  numerical¬ 
ly  integrated  to  x=0.9  and  then  the  WKB  solution  was  propa¬ 
gated  backward  to  obtain  the  comparison. 

Like  Fig  11  for  region  1,  Fig  20  is  a  plot  that  shows 
how  close  the  WKB  solutions  get  to  the  x=0.5  point  before 


Figure  20.  WKB  vs  Linear  Comparison  for  1%  Relative  Error 


Criterion  in  Region  3 


they  deviate  from  the  linear  solution  by  1%  relative  to  the 
magnitude  of  the  extreme  value  of  the  linear  solution  in 


region  3.  Fig  20  shows  that  the  best  solution  for  this 
criterion  for  e=0.1  is  one  WKB  term,  For  8=0.01,  two  terms 
came  closest  to  the  singularity.  For  8=0.001,  four  terms 
was  best . 

Fig  21  is  similar  to  Fig  12  for  region  1.  Fig  21  plots 


Figure  21.  Relative  WKB  Error  in  Last  Half  of  Region  3 


the  maximum  error  between  the  WKB  and  linear  solutions  rela¬ 
tive  to  the  magnitude  of  the  extreme  value  of  the  linear 
solution  in  region  3.  The  error  is  shown  as  a  function  of 
the  number  of  WKB  terms  used  in  the  approximation.  The 
interval  used  for  the  error  data  was  from  x=0.9  to  x=0.7  to 
try  to  eliminate  the  effect  of  the  singularity  at  X=0.5  in 
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the  WKB  approximation.  For  this  criterion,  the  best  approxi¬ 
mation  for  e=0.1  is  the  same  as  for  Fig  20  (one  term).  Now 
the  best  approximation  for  £=0.01  on  this  interval  is  five 
WKB  terms.  The  best  approximation  for  £=0.001  is  six  terms. 

Fig  22,  23,  and  24  are  examples  of  the  comparison  be¬ 
tween  the  WKB  and  linear  solutions.  Fig  22  shows  the  com¬ 
parison  for  £=0.1  and  one  WKB  term.  Fig  23  shows  the  com¬ 
parison  for  £=0.01  and  two  WKB  terms.  Fig  24  shows  the 
comparison  for  £=0.001  and  four  WKB  terms. 

The  criteria  used  in  Figs  20  and  21  are  arbitrary  and 
different  criteria  could  give  different  numbers  of  terms  as 
being  best.  Like  region  2,  care  should  be  taken  when  trying 
to  compare  one  £  to  another,  since  the  errors  in  Fig  20  and 
21  are  relative  to  different  values  for  each  £. 

Spinup  Starting  in  Region  3..  One  case  of  prolate  spin- 
up  not  examined  is  when  ji0>i2  and  13  •  This  spinup  begins  in 
region  3 .  It  was  not  examined  because  it  is  the  same  as 
oblate  spinup.  Like  oblate  spinup,  the  WKB  approximation 
will  get  worse  as  the  start  of  spinup  is  moved  closer  to  the 
WKB  singularity. 


55 


Figure  22.  WKB  and  Linear  Solutions  in  Region  3  for  e=0.1 
and  One  Term 


56 


57 


Figure  24.  WKB  and  Linear  Solutions  in  Region  3  for  e=0.001 
and  Four  Terms 
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S0  =  ±ife^2tj  <1M> 

Then  substituting  Eq  (104)  into  Eq  (76)  gives  the  same  equa¬ 
tion  as  Eq  (102)  .  Therefore,  the  WKB  solution  is  exactly 
the  same  as  the  linear  solution  for  the  axisymmetric  case. 

Also,  note  that  for  i2=i3  the  nonlinear  equations,  Eqs 
(16) -(18),  become  linear.  Now,  if  the  definition  of  k2  is 
changed  to 

^2  =  M-o  -  -*2 (105) 

the  WKB  solution  will  also  be  the  same  as  the  solution  to 
Eqs  (16) -(18)  since  X2  is  now  a  constant. 
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Conclusions 


V.  Conclusions  and  Recommendations 


In  this  thesis,  the  WKB  method  is  used  to  approximate 
the  linear  solutions  for  oblate  and  prolate  spinup  of  a 
dual-spin  spacecraft.  Based  on  the  comparisons  between  the 
WKB  and  linear  solutions,  the  following  conclusions  are 
drawn: 

1.  The  WKB  solution  is  a  good  approximation  to  the 
linear  solution  for  oblate  spinup.  For  e=0.01  and 
0.001,  the  WKB  solution  with  only  two  terms  has  a 
relative  error  of  less  than  1%  for  the  example 
problem.  For  e=0.1,  three  terms  are  required  to 
reach  this  accuracy. 

2 .  The  WKB  solution  for  prolate  spinup  is  not  prac¬ 
tical.  The  WKB  solution  contains  two  singulari¬ 
ties  around  which  the  approximation  is  not  valid. 
These  points  divide  the  WKB  solution  into  three 
parts.  To  obtain  each  succeeding  WKB  solution  the 
linear  equations  of  motion  must  be  numerically 
integrated  to  obtain  initial  conditions  for  the 
new  WKB  solution.  Also,  the  linear  solution  is 
only  valid  for  spacecraft  that  are  nearly  axisym- 
metric  and  for  very  small  perturbations  from  the 
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equilibrium  point  which  greatly  restricts  the 
conditions  under  which  the  WKB  method  can  be  used. 

3 .  The  WKB  method  provides  the  exact  solution  for  the 
axisymmetric  spacecraft.  Only  one  WKB  term  is 
required  for  the  WKB  solution  to  be  the  same  as 
either  the  solution  of  the  linear  or  nonlinear 
equations  of  motion. 

Recommendations 

The  following  are  some  recommendations  for  further 
study : 

1.  Find  some  method  to  patch  the  three  solutions  of 
the  WKB  approximation  for  prolate  spinup  together 
to  give  a  solution  that  is  not  dependent  upon  a 
numerically  integrated  solution. 

2.  Use  a  different  perturbation  method  to  approximate 
the  nonlinear  equations  of  motion.  This  would  be 
less  restrictive  than  using  the  WKB  approximations 
which  are  only  valid  where  the  linear  solution  is 
valid. 

3.  Investigate  the  optimum  range  of  e  for  the  minimum 
error  between  the  nonlinear  and  WKB  solutions  for 
spinup.  The  error  between  the  WKB  and  linear 
solution  varies  with  e  for  a  set  of  initial  condi¬ 
tions.  Generally  as  e  decreases,  the  error  de- 
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creases  depending  on  the  number  of  WKB  terms  used 
in  the  approximation.  However,  as  e  decreases, 
the  spinup  time  grows  longer  and  the  error  in¬ 
creases  between  the  linear  and  nonlinear  solution. 
This  is  shown  qualitatively  in  Fig  25. 


Figure  25.  Error  as  a  Function  of  £ 


Appendix :  Computer  Program 

Below  is  the  FORTRAN  program  that  was  used  to  generate 
the  comparison  between  the  linear  and  WKB  solutions  of  ob¬ 
late  spinup.  All  of  subroutine  HAMING  and  the  design  for 
the  main  program  and  subroutine  RHS  were  taken  from  refer¬ 
ence  7  (7:108-113) : 


program  dual spin 

common  /ham/  t,x(42,4) , £(42,4) .err (42) ,n,h, mode 
common  /dual/  12,i3,a,muO,z,zO,zpO,iterms 
double  precision  t,  x,£,err,h>12,13,e 
double  precision  z,mu0,z0,zp0 
open  (unita20,file*'C:\qpro\bln\wkb.txt' ) 

n*6 

t»0.d0 

write(*,*)  'Compare  WKB  with  exact (type t  2);  linear (type :  5).' 
read(*,*)  m 

write(*,*)  'Type  1-6  for  number  of  WKB  terms.' 

read ( * , * )  lterms 

write ( * ,  * )  'What  Is  12  ?' 

read(*, *)  12 

write(*,*)  'What  Is  13  ?' 
read(», *)  13 

write (*,*)  'What  Is  epsilon  (mu  dot)  ?' 
read ( * , * )  e 

write (*,*)  'What  is  the  initial  value  of  xl  ?' 

read ( * , * )  x(l,l) 

x(2, 1) sdsqrt (l.d0-x(l, 1) **2.d0) 

zpOeO.dO 

X(3,l)«0.d0 

write (*,*)  'What  is  the  Initial  value  of  x4  (mu)  ?' 

read ( * ,  * )  x(4,l) 

he.025d0 

write(*,*)  'This  will  produce  ' ,nlnt ( (l.d0-x(4, 1) ) / (e*h) ) , 

+  '  data  points . ' 

write (*,*)  'To  print  every  Xth  data  point,  ', 

+  'enter  I  as  an  integer  value.' 
read ( * , * )  ith 
z0«x(2,l) 
mu0>x(4,l) 
x(5,l)-x(2,l) 
x(6,l)-x(3,l) 

ii-nint ( (1 .d0-x(4, 1) ) / (e*h*dble (1th) ) ) 
nxteO 

call  haming(nxt) 
if(nxt.eq.O)  stop  99 
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do  100  1-1,11 
do  50  j -1,1th 

call  haming(nxt) 

50  continue 
call  wkb 

write (20, 1000)  t,x(m,nxt) , z, z-x(m,nxt) 
1000  format ( £8 . 3 , 2x, el6 . 9 , 2x, el6 . 9 , 2x, el£ . 9 ) 
100  continue 
c 

close (20) 
end 
c 

$ INCLUDE : ' rhsolin . for ' 

$ INCLUDE : ' bamlng . f or ' 

$ INCLUDE : 'wkbo.for' 


subroutine  rbs(nxt) 
c 

common  /bam/  t,x(42,4) ,f (42,4) , err (42) ,n,b,mode 
common  /dual/  12,i3,e,muO,z,zO,zpO,iterms 
double  precision  t,x,f,err,b,12,i3,e,mu0,z,z0,zp0 
c  These  are  the  actual  BONsi 

f (l,nxt)- (12-13) *x(2,nxt)*x(3,nxt) 
f (2,nxt)-(i3*x(l,nxt) -x(4,nxt) ) *x(3,nxt) 
f <3,nxt)--(i2*x(l,nxt) -x(4,nxt) ) *x(2,nxt) 
f (4,nxt)«e 

c  These  are  the  linearized  EOMss 
c  This  is  Z2  dot i 

f (5,nxt)«(i3-x(4»nxt))*x(6,nxt> 
c  This  is  X3  dott 

f (6,nxt)«(x(4,nxt)-i2)*x(5,nxt) 

return 

end 


subroutine  haming(nxt) 

common  /ham/  x,y(42,4) ,f (42,4) ,errest(42) ,n,h,mode 
double  precision  x,y,f , arrest, h 
double  precision  hh,xo,tol 
tol-O.OOOOOOOOOOOldO 
if(nxt)  190,10,200 
10  xo«x 

hh-h/2 . OdO 
call  rhs(l) 
do  40  1-2,4 
x-x+hh 
do  20  i-l,n 

20  y(i,l)-y(i,l-l)+hh*f (1,1-1) 
call  rhs(l) 
x-x+hh 
do  30  1-1, n 

30  y (i,l) «y(i, 1-1) +h*f (i,l) 

40  call  rhs(l) 
jsw— 10 
50  isw-1 

do  120  1-1, n 

hh«y(i,l)+h*<9.0d0*f(i,l)+19.0d0*f (i,2) 

1  -5.0d0*f (i,3)+f (i,4))/24.0d0 

if (dabs (hh-y( 1,2) ) .lt.tol)  go  to  70 

isw-0 

70  y(i,2)-hh 

hh-y (i, 1) +h* (f (1,1) +4 .0d0*f (i, 2) 

1  +  f (i,3))/3.0d0 
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if  (dabs(bb-y(i,3) ) .lt.tol)  go  to  90 
isw*0 

90  y(i,3)*bh 

bh*y(i,l)4-b*(3.0d0*f  (i,l)4-9.0d0*f  (i,2) 

1  +9 . 0d0*f (i, 3) +3 . 0d0*f (i, 4) ) /S.OdO 

If (dabs (bh-y (1,4)) .lt.tol)  go  to  110 
iSW*0 

110  y(i,4) *bb 
120  continue 
x*xo 

do  130  1*2,4 
x*x+h 

130  call  rhs(l) 

if(isw)  140,140,150 
140  jsw*jsw4-l 

if(jsw)  50,280,280 
150  x*xo 

isW*l 

jewel 

do  160  i*l,n 
160  arrest (i)*0.0d0 
nxt*l 
go  to  280 
190  jsw*2 

nxt*iabs (nxt) 

200  x-x+h 

nplamod (nxt , 4 ) +1 
go  to  (210,230) ,isw 
210  go  to  (270,270,270,220) ,nxt 
220  isw*2 

230  nnt2*mod(npl,4) +1 
nml-mod  (nm2 , 4 )  + 1 
npoemod  (nml ,  4)4-1 
do  240  1*1, n 

f (l,nm2) -y (i,npl) +4 .0d0*b* (2 .0d0*f (i,npo) 

1  -f  (l,nml)4-2.0d0*f  (i,nm2))/3.0d0 

240  y (i,npl) -f (i,nm2) -0 ,925619835d0*errest (i) 
call  rbs(npl) 
do  250  1*1, n 

y (i,npl) *(9.0d0*y(i,npo) -y(i,nm2) 

1  4-3 . 0d0*h*  (f  (i,npl)  +2 .0d0*f  (i,npo) 

2  -f (i,nml)))/8.0d0 
arrest (i) *f (i,nm2) -y(i,npl) 

250  y(l,npl) *y(i,npl)+0.0743801653d0*errest (i) 
go  to  (260,270),  jsw 
260  call  rhs(npl) 

270  nxt*npl 
280  return 
end 


subroutine  wkb 

common  /bam/  t,x(42,4) ,f (42,4) ,err(42) ,n,b,mode 
common  /dual/  i2,i3,e,mu0,x,x0,xp0,lterms 

double  precision  t,x,f,err,h 
double  precision  mu0,i2,i3,e,x,i0,sp0,k2,k3,et 
double  precision  rel,rel0,relp0,iag,iag0,iagp0 
double  precision  s0,sl,s2,s3,s4,s5 
double  precision  s00,sl0,s20,s30,s40,s50 
double  precision  S0p0,slp0,s2p0,s3p0,s4p0,s5p0 
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+  + 


c 


k2-muO-12 

k3-mu0-i3 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


sOpO-dsqrt (k2*k3 ) 

slpO- (k2-k3) / (4 .d0*k2*k3 ) 

s2p0-( (k3-k2) * (7 .d0*k2+5 .d0*k3) ) / 

+  (32.d0*(k2*k3)**2.5d0) 

s3p0- (3 .dO* (k3-k2) * (7 .dO*k2**2.dO 
+8 .d0*k2*k3+5 .d0*k3**2 .dO) > / 

(64 .dO* (k2*k3) **4.dO) 

s4p0- ( (k2-k3) * (1 .463d3*k2**3 .dO+2 . 163d3*k2**2 .d0*k3 
+  +1.989d3*k2*k3**2.d0+1.105d3*k3**3.d0))/ 

+  (2.048d3*(k2*k3)**5.5d0) 

«5p0-(3.d0*(7.07d2*k2**5.d0+5.39d2*k2**4.d0*k3 
+  +1.46d2*k2**3.d0*k3**2.d0 

+  -2.62d2*k2**2.d0*k3**3.d0 

+  -5.65d2*k2*k3**4 .dO-5 .65d2*k3**5 .dO) ) / 

+  (1.024d3*(k2*k3)**7.d0) 

do  10  1-1,2 

lf(l.eq.l)  then 

et-O.dO 

else 

et-e*t 


sO- ( (k2+k3 ) /4 .dO+et /2 .dO) *dsqrt ( (k2+et >  * (k3+et ) ) 

+  -(k3-k2) **2.d0/8 .dO*dlog (k2+k3+2 .dO*et 

+  +2.dO*dsqrt ( (k2+et) * (k3+et) ) ) 

sl-(dlog(k3+et)-dloa(k2+et) ) /4.d0 

s2— (-7.d0*k2**3.d0+27.d0*k2**2.d0*k3 
+  -9 .d0*k2*k3**2 .dO+5 .d0*k3**3 .dO 

+  +6 .d0*k2**2 .d0*et+36 .d0*k2*k3*et 

+  +6.d0*k3**2.d0*et+24.d0*k2*et**2.d0 

♦  +24 .d0*k3*et**2 .d0+16 .d0*et**3 .dO) / 

+  (48 .dO* (k3-k2 ) **2 .dO* ( <k2+et ) * (k3+et ) ) **1 . 5d0) 

s3* ( (k2-k3) * (7 .d0*k2+5 .d0*k3+12 .dO*et) ) / 

+  (64 .dO* ( (k2+et ) * (k3+et) ) **3 .dO) 

s4-(l. 105d3/ (9 .216d3* (k2-k3) **2 .dO* (k2+et) **5 .dO) 

+  +1.547d3/(9.21€d3*(k2-k3)**3.d0*(k2+et)**4.d0) 

+  +1.626d3/(1.536d4*(k2-k3)**4.d0*(k2+et)**3.d0) 

+  +1.633d3/(4.608d4*(k2-k3)**5.d0*(k2+et)**2.d0) 

+  +1 . 27d2 / ( 5 . 76d3  * ( k2 -k3 ) *  *6 . dO  * ( k2+et ) ) 

+  -1.463d3/ (9.216d3* (k2-k3) **2 .dO* (k3+et) **5 .dO) 

+  +2 . 653d3/ (9 . 216d3* (k2-k3) **3 .dO* (k3+et) **4 .dO) 

+  -4.973d3/(1.536d4*(k2-k3) **4.d0*(k3+et) **3.d0> 

+  +1.5767d4/(4.608d4*(k2-k3)**5.d0*(k3+et)**2.d0) 

+  -4.223d3/(5.76d3*(k2-k3)**6.d0*(k3+et)))* 

+  dsqrt ( (k2+et ) * (k3+et ) ) 

s5- ( (k3-k2) * (7 .07d2*k2**3 .d0+1.017d3*k2**2 .d0*k3 
+  +1 . 017d3  *k2  *k3  *  *2 . dO+5 . 6  5d2  *k3  *  *  3 . dO 

+  +3 . 192d3*k2**2 .d0*et+4 . 176d3*k2*k3*et 

♦  +2 .712d3*k3**2 .d0*et+5 .28d3*k2*et**2 .dO 
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+  +4 .8d3*k3*at**2 .dO+3 .36d3*at**3 .dO) ) / 

+  (2 . 048d3* ( (k2+at) *(k3+at) ) **6.d0) 

if(i.aq.l)  than 

sOO-sO 
BlOmal 
820-82 
830-83 
840-84 
s50«s5 
and  if 

10  continua 

if (itazms.aq.l)  than 

ralO-O.dO 

ralpO-O .dO 

imgO-sOO/a 

lmgpO-sOpO/e 

ral-O.dO 

img-sO/a 

and  if 

if (itarffls.a<z.2)  than 

ralO-slO 

ralpO-slpO 

ingO-sOO/a 

ixngpO-sOpO/a 

ral-sl 

img-sO/a 

and  if 

if (itana8.eq.3)  than 

ralO-slO 

ralpO-slpO 

img0-s00/a+a*s20 

imgpO-sOpO/e+a*82pO 

ral-sl 

lmg-80/a-fa*82 
and  if 

if (ltarms.aq.4)  than 
ral0-8l0+a**2.d0*s30 
ralp0=slp0+a**2 ,d0*s3p0 
img0-s00/e+a*s20 
imgpO-sOpO/e+e*82pO 
ral-sl+a**2 -d0*a3 
img-sO/a+a*s2 
and  if 

if (itarms.aq.5)  than 
ral0-sl0+a**2 .d0*s30 
ralpO-slpO+a**2 .d0*s3p0 
Img0«s00/a+a*820+a**3.d0*840 
imgp0-s0p0/a+a*s2p0+a**3  .d0*s4p0 
ral-sl+a**2.dO*B3 
img-s0/a+a*s2+a**3 .d0*s4 
and  if 

if (it arms. aq. 6)  than 
ral0-sl0+a**2.d0*s30+a**4.d0*s50 
ralp0-slp0+a**2 .d0*«3p0+a**4 .d0*s5p0 
img0«s00/a+a*s20+a**3 .d0*s40 
iingpO-sOpO/a+a*s2pO+8**3.dO*s4pO 
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e 


c 


c 

c 


rel»8l+a**2  .dO*B3+a**4  .d0*«5 
imgK80/a+e*B2+a**3 ,dO*B4 
•nd  if 

if (itaras.gt.6  .OR.  itarmB.lt.l)  than 

writa(*,*>  'Numbar  of  WKG  tanna  invalid I !  #tarms=' , itarms 
and  if 

z*daxp  (ral-ralO)  *  (zOMcos  (img-inaO) 

+  +(zpO-zO*ralpO) *dBin(img-imgO) /lmgpO) 

raturn 

and 
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